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ABSTRACT 

Microvariability consists in small time scale variations of low amplitude in the photometric light 
curves of quasars, and represents an important tool to investigate their inner core. Detection 
of quasar microvariations is challenging for their non-periodicity, as well as the need for high 
monitoring frequency and high signal-to-noise ratio. Statistical tests developed for the analysis 
of quasar differential light curves usually show either low power or low reliability, or both. In 
this paper we compare two statistical procedures that include several stars to perform tests with 
enhanced power and high reliability. We perform light curve simulations of variable quasars and 
non-variable stars, and analyze them with statistical procedures developed from the F-test and 
the analysis of variance. The results show a large improvement in the power of both statistical 
probes, and a larger reliability, when several stars are included in the analysis. The results from 
the simulations agree with those obtained from observations of real quasars. The high power 
and high reliability of the tests discussed in this paper improve the results that can be obtained 
from short and long time scale variability studies. These techniques are not limited to quasar 
variability; on the contrary, they can be easily implemented to other sources such as variable stars. 
Their applications to future research and to the analysis of large field photometric monitoring 
archives can reveal new variable sources. 

Subject headings: galaxies: photometry - methodsistatistical - quasars:general - techniques photometric 


1. Introduction 


Variability is an important tool to study the 
inner physics and structure of the Active Galactic 
Nuclei (AGNs). The discovery of the variable be- 
havior in the quasars ( Matthews Sz Sandaeell 1963 1 
closely followed the discovery of quasars them¬ 
selves. The time scales of the variations in the 
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optical wavelenghts range from minutes to years. 
In general, the amplitude of the variation is larger 
when observed on longer time scales and shorter 
wavelengths. Large amplitude and time scale vari¬ 
ations are relatively easy to detect, but the detec¬ 
tion becomes problematic for small amplitudes 
and time scales of minutes. Therefore, most if 
not all the microvariability studies are performed 
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We define optical microvariability as flux 
changes on time scales ranging from minutes to 
hours. The expected change in brightness for 
such events is i n the order of hundr edths of a 
magnitude (e.g. de Diego et al.l [l9^ . posing a 
challenge for detection because the amplitude of 
the variation event and the noise level are simi¬ 
lar. Moreover, quasar variability is aperiodic, and 
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therefore difficult to analyze because we cannot 
collect observations to obtain a smoothly folded 
fight curve. These limitations place a strong em¬ 
phasis on careful handling of the data and the 
use of robust and powerful statistical techniques 
for the analysis. Analyzing aperiodic variations 
has been the subject of recen t work on AGNs 
( Emmanoulopoulos et al ] l2mnl: IdePiegnl l2mfi 
2014 . the latter from now on Paper ih as well 


as yo ung stars and massive stars (iFindeisen et al 


201,511 . 


This issue was already recognized by sev¬ 
eral groups, that proposed various statistical ap¬ 
proaches to search for microvariability events. 
Thus, there are several tests that have been fre¬ 
quently used to repor t microvariatio n s. T he C- 
test as proposed b y Ijang fc Milleij ( 1997ll and 
iRomero et al.l ( 1999ll is a very simple methodology 
that gained popularity ten years ago, b ut unfor¬ 
tunately it is not a valid statistical test ( de Diegol 
2 OIOI Paper ill . The F-test is also a simple statis¬ 


tical procedure and has been replacing the C-test 
during the last years. However, the original ver¬ 
sion of the F-test requires the quasar to be com¬ 
pared with a star of the same bright ness, which is 
a condition that is difficult to meet. iHowell et al 


and Joshi et ^ ( 201lll proposed the scal¬ 
ing of the photometric errors in order to work out 
this problem. Moreover, the F-test faces another 
limitations: a lower power should be expected 
because of the violation of no rmality in data ob¬ 
tained from variable sources (IPaper Il l. The one- 
wa y analysis of va r iance (ANOVA) was proposed 
bv Ide Diego et al. ( 1998ll for microvariability de¬ 
tection, and it has been shown to be superior 
to the F-test for this purpose (|de Diegol I 2 OI 0 I 


Paper Jl . 


Another strategy that has been considered is 
the use of multiple tests (multitesting) to face 
the problem of the low reliability of some stud¬ 
ies, which can be affected by the odd behavior of 
a comparison star rather than the variability of the 
target quasar. Multitesting has been implemented 
using two diff erent tests or two different com - 
parison stars ( Joshi etldl 12011 : Hu et al. 201, ill . 
Microvariability detection is then claimed only if 
both tests, at the significance level of a = 0.01, 
agree in the rejection of the nu ll non-variability 
hypothesis. However, in lPaner^ it was shown that 
this is a low power methodology that yields unre- 
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Fig. 1.— Differential light curve examples. Open 
circles show an example of a simulation of a vari¬ 
able quasar. Filled rectangles show an example of 
a simulation of a control object. 


liable results. Nevertheless, there is the possibility 
of including several comparison stars in the light 
curve analysis in such a way that the power of the 
test would be increased, as well as the reliability. 
In this sense, the enhanced F-test, proposed in 
IPaper j involves several comparison field stars in 
a single F-test, which provides a reliable method¬ 
ology to increase the statistical power. 

In this paper we present the nested ANOVA, 
which is an improved version of ANOVA that also 
includes several stars in the quasar differential 
fight curve analysis. We compare enhanced F-test 
and nested ANOVA results using both simulations 
and examples from available real data. 

The paper is organized as follows. Section 2 de¬ 
scribes the light curve simulations, the enhanced 
F-test and the nested ANOVA methodology. In 
Section 3 we show and compare the results of 
the analysis of the simulated light curves, and 
the results obtained from real observations of the 
quasars US 3150 and IF 15498-1-203. Section 4 
presents a brief discussion and conclusions. 

2. Methods 

This study consists of one thousand light curve 
simulations for a variable quasar of magnitude 
U = 17, a non-variable control object of the same 
brightness, and 21 field bright stars between mag¬ 
nitudes 16 < U < 18, along with the comparison 
of the simulations with real observations. Figure[T] 
shows an example of the simulations of differen- 
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tial light curves: a variable quasar and a non¬ 
variable control star with the same photometric 
error. Each differential light curve consisted of 
35 observations obtained by subtracting the raw 
magnitudes of a reference star from the target to 
eliminate the changing-conditions effects between 
exposures. This reference star is the brightest, 
non-saturated star in the held. 

In contrast, the ANOVA procedures use only 
reference stars, but the test is performed between 
different groups of differential photometry obser¬ 
vations of the same target rather than between 
th e target and a comparison star light curves. As 
in lPaperl . in the case of the variable quasars, the 
light curves were the result of a random walk func¬ 
tion with steps normally distributed with mean 0 
and standard deviation 0.006 [A/’(0,0.006) RW], 
and a normally distributed photometric error of 
e ~ 0.007. The photometric errors for the set of 
21 held stars are also normally distributed. All 
the photometric errors correspond to 60s expo¬ 
sures in the V band with the Harold L. Johnson 
1.5 m telescope at the Observatorio Astronomico 
Nacional in Sierra San Pedro Martir (Mexico), and 
they were calculated using the simulator for this 
instrument 0 

As detectors become larger and detector ar¬ 
rays more popular, there is a tendency towards 
larger held sizes and thus images that include more 
bright stars. This situation makes it necessary to 
perform realistic simulations. To deal with this 
necessity, the magnitudes for the 21 held stars 
were obtained using the star distribu tion in the V 
band for a galactic latitude = 60° ( Allen fc CoxI 
2000L Table 19.11, pp. 482-3). The star distri¬ 
bution in magnitudes was approximated using a 
second order polynomial: 

logNv = -4.026 -f 0.56071/ - 0.00984E^ (1) 

where Ny is the cumulative function of the num¬ 
ber of stars brighter than magnitude V in one 
square degree of the sky. The V magnitudes for 
the 21 reference and comparison stars were ob¬ 
tained applying the inversion method of the cu¬ 
mulative function to the uniform distribution on 
the Ny values in the range 16 < 1/ < 18, i.e. up 
to one magnitude of difference with respect to the 


Table 1: Star parameters for simulations. 


Index 

V 

Error 

Index 

V 

Error 

1 

17.727 

0.011 

11 

16.119 

0.004 

2 

17.998 

0.013 

12 

16.655 

0.006 

3 

17.798 

0.012 

13 

17.392 

0.009 

4 

17.173 

0.008 

14 

17.336 

0.009 

5 

17.351 

0.009 

15 

17.850 

0.012 

6 

16.709 

0.006 

16 

17.493 

0.010 

7 

17.761 

0.011 

17 

17.856 

0.012 

8 

16.633 

0.006 

18 

17.077 

0.007 

9 

16.185 

0.004 

19 

17.703 

0.011 

10 

16.942 

0.007 

20 

16.871 

0.007 




21 

17.064 

0.007 


quasar. This range of magnitudes ensures that the 
images of the stars are suitable for differential pho¬ 
tometry with the quasar, i.e. images are not suf¬ 
fering either overexposure or underexposure. The 
parameters for the stars, in the order in which they 
appear in the simulations, are shown in Table [T] 

2.1. Enhanced F-test 

The simulations were an alyzed us ing the en¬ 
hanced F-test, described in Paper 1 Until now, 
the F-test has been implemented with a single 
comparison star, whose differential light curve is 
compared with the differential light curve obtained 
from a target object using the ratio of their re¬ 
spective variances. It is important that the mean 
brightness of both the comparison star and the 
target object are matched to ensure that the pho¬ 
tometric errors are equal, or at least that these 
errors are corrected to accoui it for the differences 
in br ightness, as propos ed bv iHowell et al. ( 1988h 


and iJoshi et al.l ( 2011 ). The enhanced F-test 


^Developed by Alan Watson, V. 1.2.0: 
http: //www. astrossp. unam. nix/'/,7Eresast/wat sonBCh/\f otometria 


makes use of several comparisons stars, and it con¬ 
sists in transforming the comparison star differ¬ 
ential light curves to have the same photometric 
noise as if their magnitudes matched exactly the 
mean magnitude of the quasar under study. Note 
that the enhanced F-test needs to use a reference 
star in order to obtain the differential light curves. 
Thus the number of available bright stars to per¬ 
form the test is reduced by one, and the test is 
performed using only the rest of bright compari¬ 
son stars. 

-ccd.html 
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Fig. 2.— Enhanced F-test steps. Panel (a) shows the magnitudes of the 21 simulated stars and their 
associate errors obtained with the San Pedro Martir instrument simulator. Panel (b) shows the differential 
light curve means and their associated errors for 20 comparison stars using the brightest star as reference; 
the dashed line indicates the exponential fit to the errors. Panel (c) shows the differential light curves for 
the 20 stars combined in a single set of 20 x 35 = 700 indexed observations (the star bins, of 35 observations 
each, are identified in the upper axis);after subtracting each star mean and transforming the errors to the 
same level as the quasar (in our case using the exponential fitting function), the 20 differential transformed 
light curves are stacked (lower red light curve). Panel (d) zooms on the stacked light curve shown in the 
previous panel. 


Figure [2] helps us to understand the algorithm 
to perform the enhanced F-test both in our sim¬ 
ulations and in a real case. Figure [5^ shows the 
error and the V magnitudes of the stars in a given 
field; the errors were obtained from the San Pe¬ 
dro Martir instrument simulator for the 1.5 m tele¬ 
scope. As in a real situation, we only need instru¬ 
mental magnitudes to build the differential light 
curves. We use the brightest star as our reference 


star to minimize errors; in our simulations this star 
is labeled with the number 11 (Tabled]). Figure[2)D 
shows the errors vs. differential magnitudes. We 
obtained these errors and magnitudes empirically 
from the simulated differential light curves, using 
their standard deviations and means, respectively. 
The dashed line shows a fit to the differential mag¬ 
nitudes and errors data. In this paper we have 
used an exponential fit (see Section n. It is not 
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important how we describe the dependence be¬ 
tween instrumental erro rs and magnitudes , or even 
if we us e an empirical d.Toshi et al. 2011 1 or the¬ 


oretical (|Howell et alJ 1198811 relationship, as long 
as the fits are reasonable. The point that matters 
is that we can estimate the differential photomet¬ 
ric error for our target quasar. For a quasar of 
magnitude Vq = 17 and the bright reference star 
11 (Vii = 16.119), the differential magnitude is 
AVq =Vq — Vii = 0.881, and the exponential fit 
yields an error cav = 0.008. Now we can trans¬ 
form the comparison star differential light curves 
to test the quasar variations. Figure [2j: shows the 
differential light curves for each comparison star. 
We subtract the mean for the z-star from the re¬ 
spective light curve, and multiply the resulting 
dataset by the ratio of the quasar and the star 
fitted errors: eAv{Q)/ s-av^SU). Then we have 
transformed the original comparison star differen¬ 
tial light curve to values that can be compared 
with the quasar data. Stacking all the transformed 
comparison star light curves, we obtain a single 
dataset to perform the enhanced F-test (longer 
light curve at the bottom of the plot). Note that 
in the F-test the variances, rather than the means, 
are compared; our stacked dataset has zero mean, 
but it does not affect our results (an arbitrary con¬ 
stant value, for example the quasar mean differen¬ 
tial magnitude, may be added to this dataset for 
graphical or other cosmetic purposes). Figure [2ji 
shows the stacked dataset with more detail. 


2.2. Nested ANOVA 

We also analyzed the simulations using the 
nested ANOVA probe. The tests of the ANOVA 
family require groups of replicated observations. 
In our case, these groups are arranged by five ob¬ 
servations of a given quasar in a time lapse of 
around 5 to 10 minutes for telescopes with aper¬ 
tures larger than 1 m, so that the quasar bright¬ 
ness can be considered constant for the observa¬ 
tions collected in the same group. This is the 
same procedure used in pr evious real observa¬ 


tion ( de Diego et al. 1998 : Bamfrez et al. 2004 


2 OO 9 II and simulations ( de DiegcJ 2010l Paper jb 


Of course, microvariations shorter than the lapse 
time within the group observations would not be 
detected by ANOVA, but such phenomena, known 
as spike s, has been seldom reported in the lit- 
erature ( Sagar et al. 1996[ de Diego et al.l 1998; 


Krishna et al. 2000t Stalin et al. 2004ll . 


The nested ANOVA, as implemented in this 
paper, consists in the study of the variances at 
three different experimental levels. The first level 
is the one in which we are really interested: the 
differences between the groups, that are a signa¬ 
ture of the quasar variability. The second level 
corresponds to the differences between the obser¬ 
vations due principally to shot noise and sky sub¬ 
traction. The third level accounts for the variance 
between the tests caused by the different refer¬ 
ence stars. A complete description of the nested 
ANOVA procedure is given in Appendix [A] Cal¬ 
culations have been performed using the R lan¬ 
guage, and the nested ANOVA probes were car¬ 
ried out using the aov function included in the 
stats package. The complete code used for the 
nested ANOVA analysis applied to real data can 
be found in Appendix IbI 


Note that ANOVA family tests compare the dis¬ 
persion of the individual differential magnitudes of 
the quasar within the groups, and the dispersion 
between the groups, without using any compari¬ 
son star. Instead, all the bright stars are utilized 
as reference stars to build distinct differential light 
curves, and thus we will always be able to include 
one more star in the nested ANOVA analysis than 
in the enhanced F-test. 


3. Results 


3.1. Simulations 


The results in Paper j showed that ANOVA 


has more power than the F-test to detect vari¬ 
ability when only one comparison star is taken 
into account in the analysis. Moreover, the power 
of the F-test was reduced by the probably non¬ 
normal distribution of the variable quasar data. 
This power loss can be compensated by using the 
enhanced F-test, that includes more comparison 
stars in the analysis. Here we present the re¬ 
sults of increasing the number of stars using nested 
ANOVA that lead to a test power improvement, 
and we compare the outcomes obtained using both 
the nested ANOVA and the enhanced F tests. 


The stars in the simulations have different mag¬ 
nitudes and errors, as explained in Section!^ Fig- 
ure[3]shows the number of false detections or Type 
I errors, and the number of detections. The stars 
are included in the same order as in Table H] to 
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Number of Bright Stars Number of Bright Stars 

Fig. 3.— Test outcomes for 1000 light curve simulations. Results for the enhanced F-test and nested ANOVA 
are shown as dashed line connected squares, and solid line connected circles, respectively. Tests have been 
performed at the signihcance levels a = 0.01 [panels (a) and (b)] and a = 0.001 [panels (c) and (d)]. Type 1 
errors obtained from the control object are shown in panels (a) and (c), and the number of null hypothesis 
rejections for the varying quasar are shown in panels (b) and (d). The magnitudes and photometric errors 
for the stars have been obtained from a realistic star held at galactic latitude = 60°, as explained in the 
text. 

reproduce the characteristics of the available ob¬ 
jects as the held size increases (i.e., the less nu¬ 
merous bright stars will show up as the held size 
increases). However, note that the hnal result for 
the test with a given number of stars is indepen¬ 
dent of the order in which these stars are intro¬ 
duced for both the enhanced F-test and nested 
ANOVA. The reason why the number of stars are 
different for both tests is that both the F-test 
and the enhanced F-test require an extra refer¬ 
ence star for differential photometry and the other 
stars for comparison, while for the nested ANOVA 
test all the stars are used as reference (no compar¬ 
ison stars are necessary). As a result, given a set 
of bright stars in the quasar held, we will always 
have one more star to perform the nested ANOVA 
probe than the enhanced F-test. To account for 


this effect, the niche for the bright star number 
1 in Figure [3] is empty for the enhanced F-test, 
but not for the nested ANOVA. The number of 
bright stars is directly proportional to the size of 
the held. For a galactic latitude = 60°, we ex¬ 
pect around 0.13 stars per square arcminute, i.e. 
around 2 stars in a 4' x 4' detector. 

As the number of stars included in the analy¬ 
sis increases, the tests are more and more pow¬ 
erful tending to an upper asymptotic value. We 
see in Figure [3] that in our simulations the num¬ 
ber of detections for the enhanced F-test is larger 
than the number for the nested ANOVA for a num¬ 
ber of eight or more stars. The enhanced F-test 
starts with very few detections for a few and dim 
stars (that dominate the star distribution), but the 
number of detections increases dramatically and 
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overcomes the detection for the nested ANOVA as 
the number of stars increases and bright stars are 
introduced into the analysis. 

The behavior of both the nested ANOVA and 
the enhanced F-test can be explained as the re¬ 
sult of the test designs and the data attributes. In 
the case of nested ANOVA, all the stars are used 
as reference stars, and all of them have the same 
weight in this probe. However, the results of the 
enhanced F-test depend very much on the photo¬ 
metric errors of the reference star. These errors 
affect the precision of the whole differential pho¬ 
tometry. As the dimm stars dominates the star 
distribution, we need large fields to include bright 
stars. We can see this effect on the enhanced F- 
test detection curves in Figures [Sjj and|3ji. When 
a bright star such as the star number 6 becomes 
the reference star, the number of the enhanced 
F-test and nested ANOVA detections are compa¬ 
rable, and when the even brighter stars 8 and 9 
successively became the reference stars (and star 
6 becomes a comparison star), the enhanced F- 
test overcomes nested ANOVA. 


The inclusion of comparison stars brighter and 
dimmer than the quasar allows a better fit for the 
photometric errors as a function of the differential 
magnitudes. If all the stars are dimmer than the 
quasar, the exponential fit that we have used must 
be extrapolated to the quasar brightness, which is 
less accurate than interpolating the error. In fact, 
we have found that this inaccuracy produces too 
many type I errors. This may be an inconvenience 
in the cases of few field stars or of bright quasars. 
If comparison stars brighter than the quasar are 
not available, at least for simulated data it is safer 
to give more weight to the star with the bright¬ 
ness which is more similar to that of the quasar, 
although this procedure reduces the power of the 
enhanced F-test. In a real situation, the errors 
may be indi vidually estimated using another pro - 
cedure le.g. iHowell et al.l 1988t Joshi et al. 2011 ). 
but generalizing such a procedure for simulated 
data is beyond the scope of this paper. 

As the tests with more stars are not fully inde¬ 
pendent of the test with less stars (for example, 14 
of the 15 stars included in the 15 star tests are the 
same as in the tests for 14 stars), the Type I errors 
are not random, but show a memory of the errors 
obtained in the previous tests. Finally, the effect 
size estimate (for example, the F value of the ra- 



Comparison Star 


Fig. 4.— Two examples of the enhanced F-test 
output. Panel (a) shows the F value calculated as 
the variance ratio between the simulated quasar 
light curve and the stacked star light curves for a 
clearly variable quasar (filled squares) and a less 
variable one (hollow circles). Panel (b) shows the 
logarithm of the probabilities that the data from 
the earlier quasar have been drawn from a non¬ 
variable object. The probabilities show a decreas¬ 
ing trend, as expected from the power increase 
when more stars are added in the analysis. Panel 
(c) is like the panel (b) but for the later quasar. 
Although the F values show very small variations 
for this object, and the probabilities show the ex¬ 
pected trend as in panel (b), the tests are just 
in the limit of detection at the 0.001 level (dashed 
line), and tend to fluctuate between detections and 
non-detections. 
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Fig. 5.— Monitored light curves for quasars US 3150 (November 13, 1996), and IE 15498+203 (May 12, 
1997) and a comparison star. Panels (a) and (b) show the US 3150 raw and binned differential light curves 
(upper points) and a comparison star (lower points). Error bars in panels (b) and (d) correspond to the 
standard error measured within each group of observations. 


tio of variances in the enhanced E-test) is fairly 
constant when adding more stars in the analysis. 
Although F approaches an asymptotic value and 
the probabilities decline smoothly with the num¬ 
ber of stars (Figure|4l3), for a given light curve the 
analysis may yield fluctuations in the probabilities 
obtained with different number of stars. In those 
cases when the F value is close to the critical limit 
of detection, the tests may fluctuate between de¬ 
tections and non-detections as shown in FigurelH:. 


3.2. Observations 

Data discussed in this section correspond to the 
radioquiet quasars US 3150 and IE 15498+203, 
and h ave been already reported in de Diego et al.l 
( 1998l l. Here we present a new analysis including 
more field stars to compare the results obtained 
with nested ANOVA and the enhanced E-test. 


The quasar light curves are shown in Figure [S] 
Observations of US 3150 [V ~ 16.8) were obtained 
during about 6h of monitoring on November 13, 
1996 with the //13.05 1.5 m telescope at San Pe¬ 
dro Martir Observatory (Baja California, Mexico) 
equipped with a Thomson THX 31156 CCD of 
1024 X 1024 pixels of 19 x 19 /im and Metachrome 
II coated. Quasar US 3150 showed clear microvari- 
ability at a level of signiflcance of a = 0.001, an 
amplitude of Am ~ 0.1 mag. The data consist 
of 7 groups of 5 observations with exposures of 
60s in the V band. There were five stars in the 
3.'4 X 3f4 held suitable for our purpose, but at 
the time when we performed the observations, we 
had no plans to use all these stars in our analy¬ 
sis. Hence, the brightest star (more than 3 magni¬ 
tudes brighter than the quasar) was overexposed 
in the two last groups of observations. Therefore, 












Table 2: Test results for US 3150 data. 


N. Stars 


US 3150 



Control 


E-test 

ANOVA 

E-test 

ANOVA 

E 

P 

E 

P 

E 

P 

E 

P 

1 

2.240 

2.2E-2 

13.44 

1.8E-5 

0.622 

0.87 

0.995 

0.43 

2 

2.573 

2.9E-3 

15.71 

5.8E-6 

0.671 

0.85 

1.244 

0.33 

3 

2.450 

2.0E-3 

17.52 

2.6E-6 

0.639 

0.89 

1.123 

0.37 

4 

2.525 

8.IE-4 

15.34 

6.9E-6 

- 

- 

1.767 

0.18 

5 

- 

- 

17.06 

3.1E-6 

- 

- 

- 

- 


Note. —The number of degrees of freedom for the enhanced T’-test is 24 for the numerator and 22, 46, 70 and 94 for the 
denominator, for 1, 2, 3, and 4 stars respectively. For the nested ANOVA, the number of degrees of freedom is always 4 for the 
numerator and 20 for the denominator. 


only for demonstration purposes, we have dropped 
these two groups and we performed the analysis 
using the first 5 sets of observations, correspond¬ 
ing to about 4h of monitoring. 

Observations of IE 15498-1203 (V ~ 17.12) 
were obtained during 4.5 h of monitoring on 
May 12, 1997 with the same telescope, but this 
time equipped with a Tektronix TK1024AB of 
1024 X 1024 pixels of 24 x 24 ^m, thinned and 
Metachrome II coated. Quasar IE 15498-1-203 
showed only weak evidence of variability at a = 
0.01, an amplitude of Am ~ 0.06 mag. The data 
consists of 5 groups of 5 observations with expo¬ 
sures of 60s in the V band. There were four stars 
in the 413 x 413 field suitable for our purpose, but 
photometry was affected by the scattered light 
from the bright star HD 142109 {V = 8.94) at 
around 1' from the quasar, which spoils the ob¬ 
servations of some stars in the third group. Con¬ 
sequently, we have dropped this group from our 
analysis. 

Tables! and [3] show the results obtained with 
the enhanced E-test and nested ANOVA for 
US 3150 and IE 15498-1-203, respectively. In both 
tables, column 1 indicates the number of stars 
used in the analysis. Columns 2 and 3 show the 
F value and the associated probability for the en¬ 
hanced E-test. Columns 4 and 5 indicate the same 
values for the nested ANOVA probe. Columns 6 
and 7, and 8 and 9, show the same values for the 
control star. The number of stars in the enhanced 
E tests is one less than in the corresponding nested 


ANOVA probes because the brightest star in the 
field was used as reference. The number of stars 
in the control tests is one less than in the target 
quasar because the star most similar in brightness 
to the quasar was used as control. 

3.2.1. Enhanced F-test 

We have used the brightest star as reference to 
build the photometric differential light curves. To 
perform the enhanced E-test it is necessary that 
the errors of the comparison stars are transformed 
to the same level as the errors of the quasar. In 
principle, any function that reasonably fits the dif¬ 
ferential light curve standard deviations and the 
mean magnitudes relationship can be used. How¬ 
ever, in another work (Polednikova et ah, in prepa¬ 
ration) we have fitted exponential curves to a large 
set of comparison stars in larger CCD fields, and 
find that such a function accurately describes this 
relationship. The function is of the form: 

Si = So -l- Ae™% (2) 

where Si is the light curve standard deviation and 
rui is the light curve mean magnitude for the ith 
star, while sq and A are the parameters to be fit¬ 
ted. Notice that fitting two parameters to the 
comparison stars data reduces in two units the 
number of degrees of freedom for the comparison 
stars of the enhanced E-test. Figure |6] shows the 
relationship between differential light curves stan¬ 
dard deviation for the comparison stars and the 
quasars (excluded from the fit), and their mean 
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differential magnitudes. The exponential curve fits 
the data for the comparison stars up to a precision 
of 0.001 mag. 

The data for the quasars shown in Figure [5] are 
clearly far away from stars of comparable bright¬ 
ness and the fitted exponential curve, indicating 
that their standard deviation may be dominated 
by variability rather than photometric errors. Yet, 
the F'-test with a single comparison star fails to 
detect variations in US 3150 at the level of sig¬ 
nificance a = 0.01, and 4 comparison stars are 
needed to reach a = 0.001 (Table[2]). In the case of 
IE 15498-1-203, the enhanced F-test only reaches 
the level of significance a = 0.01 when using three 
stars. The F statistics for the enhanced F-test 
on both quasars shows small fluctuations but a 
monotonic decrease in the p value. Finally, the 
enhanced F-test for the control stars also shows 
small fluctuations (in US 3150) for both the F 
statistics and the p values, but in no case do the 
tests show any evidence of variation. 


3.2.2. Nested AN OVA 


Results of the nested ANOVA probes presented 
in Tables [2] and [3] show fluctuations for the F 
statistics and the p values for both the quasars 
and the control stars. In the case of US 3150, all 
the tests yield results below the significance level 
a = 0.001, but for IE 15498-1-203 the p values are 
only below a = 0.01. In both cases, the control 
stars do not show any evidence of variations. 


The results obtained by both the enhanced 
F-test and nested ANOVA agree that US 3150 
shows strong evidence of microvariability, while for 
IE 15498- 1-203 the evidence is no t conclusive, as re¬ 
ported in de Diego et al.l ( 1998ll . They also agree 
with the results from the simulations, in the sense 
that nested ANOVA has more power to detect mi¬ 
crovariations than the enhanced F-test when there 
is a limited number of available field stars. 


4. Conclusions 

In this paper we have presented the nested 
ANOVA test and have compared the power of 
this pro be with the enhanced F-test presented in 
IPaner j . These tests use several field stars to an¬ 
alyze the light curve of a target quasar. Both 
the enhanced F-test and nested ANOVA are very 
powerful to detect photometric variations in differ- 




Fig. 6.— Light curve standard deviations vs mean 
differential magnitudes dependence. The standard 
deviation of the light curves yield error estimates. 
Panel (a) shows four comparison stars (squares) 
in the quasar US 3150 (circle) field that have been 
used to ht an exponential curve to the relationship 
between the standard deviation and the mean dif¬ 
ferential magnitude (dashed line). Panel (b) is like 
the previous panel but for three comparison stars 
(squares) in the field of the quasar IE 15498-1-203 
(circle). These figures show that the data disper¬ 
sion for the quasar was much larger than the dis¬ 
persion for their respective comparison stars, in¬ 
dicating possible variations. 


ential light curves, and they asymptotically tend 
to larger power values than their respective single 
star counterparts (the F-test and ANOVA). The 
nested ANOVA probe shows a superior perfor¬ 
mance when the number of available bright stars is 
limited, but the enhanced F-test surpasses nested 
ANOVA when there are several bright stars, some 
of them brighter than the target object. This is 
a consequence of using a bright reference star by 
the enhanced F-test that diminishes the errors of 
the differential photometry. 
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Table 3: Test results for IE 15498+203 data. 


N. Stars 


IE 15498+203 



Control 


F-test 

ANOVA 

F-test 

ANOVA 

F 

P 

F 

P 

F 

P 

F 

P 

1 

3.098 

1.2E-2 

5.231 

l.OE-2 

1.252 

0.32 

0.361 

0.78 

2 

2.384 

1.2E-2 

7.198 

2.8E-3 

0.963 

0.52 

0.467 

0.71 

3 

2.477 

4.6E-3 

5.808 

7.0E-3 

- 

- 

1.867 

0.18 

4 

- 

- 

6.829 

3.6E-3 

- 

- 

- 

- 


Note. —The number of degrees of freedom for the enhanced F-test is 19 for the numerator and 17, 36, and 55 for the 
denominator, for 1, 2, and 3 stars respectively. For the nested ANOVA, the number of degrees of freedom is always 3 for the 
numerator and 16 for the denominator. 


In our simulations, the power increases with the 
number of stars included in the analyses from ap- 
proximately 40% to 65% for the nested ANOVA 
probes, and from 5% to 90% for the enhanced 
F-test, both at the significance level of a = 
0.001. These results agree with those obtained 
from the observed light curve of quasars US 3150 
and IE 15498+203, for which the probabilities 
tend to lower values as the number of stars in¬ 
creases, while nested ANOVA yields lower proba¬ 
bilities than the enhanced F-test due to the lim¬ 
ited number of field stars. In this paper, we have 
demonstrated that the power of microvariability 
studies increases significantly by including several 
stars in the analysis, and that this procedure is 
an important improvement with respect to all the 
other tests that have been used in the quasar mi¬ 
crovariability literature, from the C-test to the sin¬ 
gle star F-test and ANOVA, as well as multitesting 
techniques. 

From these results, we conclude that both the 
enhanced F-test and the nested ANOVA probe are 
robust techniques to detect quasar microvariabil¬ 
ity. They do not only add power to the microvari¬ 
ability detection, but by introducing several stars 
in the light curve analysis, ensure more reliable re¬ 
sults than those obtained when using a single star. 
In small fields with a limited number of stars, we 
encourage the use of the nested ANOVA rather 
than the enhanced F-test because under these cir¬ 
cumstances nested ANOVA is more powerful, and 
it can include one more bright star in the analysis. 

The use of both the enhanced F-test and nested 


ANOVA in future research will improve the de¬ 
tections of microvariability in quasars. This im¬ 
provement will eventually yield a better under¬ 
standing of the physics involved in microvariability 
phenomena and the central quasar’s engine. The 
enhanced F-test can be easily applicable to select 
varying sources candidates (i.e. quasars and vari¬ 
able stars) in regions visited several times during 
large photometric surveys, such as the the Legacy 
Survey Stripe 82 from the Sloan Digital Sky Sur¬ 
vey. The wide field images in such surveys, that 
include a large number of stars, allow to build a 
mean reference star with sources different to the 
stars used for comparison in the enhanced F-test, 
thus ensuring the independence of the reference 
mean star and the comparison stars. 
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Fig. 7.— nested ANOVA diagram. The goal is to test differences between groups of observations. Each 
group consists of 5 images that include the quasar and a number of stars. We can obtain a differential 
photometry measurement of the quasar from each star. 


A. Nested ANOVA 


The nested ANOVA analysis presented in this paper involves simulations of a = 7 groups, each of them 
consisting of 6 = 5 images (35 in total) of a target quasar and a number n < 21 of held stars that can be used 
for replicate measurements of the quasar differential photometry. In practice, the images contained in a given 
group were obtained in a time lapse that is short (less than 10 min) in comparison with the microvariability 
time scale of the quasar, while the gap between groups of observations that will be compared is arbitrary. 


This methodology implies n pseudoreplicates per image to yield a total of 35n readings of the quasar 
differential photometry. We use the term pseudoreplicates for the n readings of a given image to note that 
these readings are not independent because we use the same quasar observation (the same image) to estimate 
the quasar differential photometry with n stars. Images are the experimental units of our study; as we have 
35 images (34 degrees of freedom in total), no matter how many stars we use in our analysis, the degrees 
of freedom are divided into 6 for the (7) groups of images and 28 (= 34 — 6) for errors. The number of 
replicate measurements that we perform on a given image will improve the precision on that image, and the 
variability among the images from a given group yields a measurement of the variability within that group. 

Figure [7] shows the schematic representation of the observational methodology: an arbitrary number of 
groups of observations (in our example 7), 5 images of the quasar in each group, and n differential photometry 
measurements from each image, using n differe nt reference stars. In fact, this is a two-stage nested design, 
probably the simplest nested ANOVA. Following lMontgomerv ( 20I3L chapter 14), the linear statistical model 
is expressed by: 

Uijk — M 4” T ^j{i) 4” (-^4) 


where yijk corresponds to the measurement of the quasar differential photometry using the reference star 
k = 1, 2,..., n, located in the image j = 1, 2,..., 6 of the i = I, 2,..., a group of observations. The true 
mean of the quasar light curve is denoted by p, and the deviation from the true mean of the fth group is 
denoted by 7 ^. The deviation of the jth image of the quasar with respect to the mean of the fth group is 
denoted by a;j(q, where the subscript j(i) indicates that the j observation is nested under the ith. group. 
Finally, Sijk corresponds to the error term. 

The model described by equation lAII yields a sum of (deviation) squares that, after some manipulation 
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( Montgomery 2013h . can be expressed as: 


a b n a a b a b n 

- yf = bnY^{y, - yf + nYY^y^i ~ ^YYY^y^i^ ~ ^ 

2=1 j — 1 k—1 i—1 i—1 j=l i—1 j—1 k—1 

where the horizontal line over a letter indicates the mean value [e.g., yi = {Y^^j=i^k=iyijk)/bn]. Equa¬ 
tion (IA2l can also be expressed symbolically as: 


SSt = SSg + SSoiG) + SSe, ( A3) 

where SSt denotes the total sum of squares, SSg the sum of squares due to groups, S'S'o(G) the sum of 
squares due to the nested observations in groups, and SSe the residual sum of squares due to errors. The 
degrees of freedom for each square sum is: abn — 1 for SSt, a — 1 for SSg, ^{b — 1) for SSo{g) , s-nd ab(n — 1) 
for SSe- Dividing a sum of squares S'S' by the respective degrees of freedom v yields the mean square 
(MS = SS/v). The ratios between mean squares are distributed as F. 

In our simulations, the groups of observations are drawn from an infinite pool of instants in the continuous 
light curve of a target object, and the observations in each group are drawn from an infinite pool of possible 
images, thus both groups and observations constitute random effects. In this case we assume that is 
NID(0,cr^) and is NID(0,tT2.); where NID(0,cr^) means normally and independently distributed with 
mean 0 and variance Then we test Mg : = 0 by MSg/MSo{g) and Mg : ctJ = 0 by MSo{g)/MSe- 

The residuals for a two-stage nested design are given by 

^ijk — yijk yij, ( A4) 

where yij are the individual image averages, and the fitted value is: 

Vijk =y + li+ i^j(i) (A5) 

For random effects, the analysis of variance method allows to estimate the variance components cr|, 
and cr^ using the expected mean squares: 


= MSe 


(A6) 


^2 MSo{g) — MSe 

(7 a — - 




MSg — MSo(g) 

bn 


(A7) 

(AS) 


B. Nested Analysis of Variance R code 

# Nested ANOVA 

# (c) J. A. de Diego - November 2014 

# Input file example (file without headers) 


# TIME 

QUASAR 

SI 

S2 

S3 

S4 

S5 

GROUP 

# 03:08:50 

21.504 

17.389 

17.851 

18.911 

21.187 

21.577 

1 

# 03:10:44 

21.557 

17.347 

17.882 

18.951 

21.232 

21.655 

1 

# 03:12:34 

21.488 

17.400 

17.840 

18.913 

21.207 

21.597 

1 

# 03:14:24 

21.481 

17.319 

17.844 

18.908 

21.178 

21.593 

1 


14 









# 03:15:54 

21.484 

17.348 

17.841 

18.919 

21.176 

21.577 

1 

# 04:31:24 

21.601 

17.354 

17.877 

18.938 

21.232 

21.611 

2 

# 04:33:29 

21.552 

17.338 

17.872 

18.930 

21.222 

21.597 

2 

# 04:34:54 

21.635 

17.368 

17.895 

18.966 

21.217 

21.647 

2 

# 04:36:19 

21.583 

17.334 

17.863 

18.925 

21.233 

21.607 

2 

# 04:37:44 

21.581 

17.363 

17.894 

18.963 

21.240 

21.651 

2 


# The number of columns for the stars is arbitrary 

# Magnitudes for the quasar and stars are as generated byt IRAF 

# routines, i.e. unprocessed (not differential or extinction corrected) 

# Columns TIME and GROUP are not used in this code, but probably needed 

# for figures and other computations. 

# This code assumes that the last star corresponds to the control. 

# Output 

# nest.anova.sum.qv Contains a list with the results of tests for the quasar 

# nest.anova.sum.qn Contains a list with the results of tests for the control star 

# prob.anova.qv Is a vector with the probabilities of tests for the quasar 

# prob.anova.qn Is a vector with the probabilities of tests for the control star 


# Read data from "mydata.dat" and prepare variables 
(test <- read.table("US3150_V3.dat")) 

s <- length(test) -2 # Number of sources in the field (including the quasar) 

n <- dim(test) [1] # Data points in the light curve 

ng <- 5 # Number of points in each group 

g <- n/ng # Number of groups 

qlc <- test 

colnaunes(qlc) <- cC'Time","q",pasteO("s",seq(l:(s-1))),"Group") 
dqlc <- data.frame( 

qv = as.vector(as.matrix(qlc$q-qlc[,c(pasteO("s",seqd:(s-1))))])), # Variable quasar 

qn = as.vector(as.matrix(qlc$s5-qlc[,c(pasteO("s",seqd:(s-1))))])), # Control star 

St = rep(l:(s-l), each=n), # Star factor 

ag = rep(l:g, each=ng), # ANOVA group factor (or read column Group) 

ao = rep(l:ng) # ANOVA group element factor 

) 

# Allocate matrices 

nest.anova.sum.qv <-vector("list") 
nest.anova.sum.qn <-vector("list") 
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prob.anova.qv <- matrix(,nrow = 1, ncol = (s-1)) 

colnames(prob.anova.qv) <- paste(’cs1:(s-1), sep="") 
prob.anova.qn <- matrix(,nrow = 1, ncol = (s-1)) 

colnames(prob.anova.qn) <- paste(’cs1:(s-1), sep="") 

# Perform test 
for(j in 2:(s)) { 

# Temporal quasar variables 
dlcj <- dqlc[dqlc$st < j,] 

ym <- as.vector(tapply(dlcj$qv, 

list(dlcj$ag, dlcj$ao), 

mean)) # vector of observation (image) means 
tm <- factor(as.vector(tapply(as.numeric(dlcj$ag), 

list(dlcj$ag, dlcj$ao), 
mean))) # group factors 

# Test and probabilities for the quasar 

nest.anova.sum.qv[[(j-1)]] <- summary(aov(ym~tm)) 

prob.anova.qv[j-l] <- nest.anova.sum.qv[(j-1)] [[!]] [[!]] [,"Pr(>F)"][[!]] 

# Temporal control variables 
ym <- as.vector(tapply(dlcj$qn, 

list(dlcj$ag, dlcj$ao), 

mean)) # vector of obsevation (image) means 
tm <- factor(as.vector(tapply(as.numeric(dlcj$ag), 

list(dlcj$ag, dlcj$ao), 
mean))) # group factors 

# Test and probabilities for the quasar 

nest.anova.sum.qn[[(j-1)]] <- summary(aov(ym~tm)) 

prob.anova.qn[j-1] <- nest.anova.sum.qn[(j-1)][[!]] [[!]] [,"Pr(>F)"][[!]] 

} 

# Remove last element for the control star results (it is a test with itself) 
nest.anova.sum.qn[[(j-1)]] <- NULL 

prob.anova.qn <- prob.anova.qn[-(j-l)] 
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